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ABSTRACT 


The  early  body-wave  coda  of  far-regional  events  (14°-  29  degrees)  contains  triplicated  arrivals  from  upper-mantle 
discontinuities  that,  when  properly  identified,  ean  improve  seismic  monitoring  functions.  However,  far-regional 
seismograms  are  typically  under-utilized  because  along-path  heterogeneity  and  phase  interactions  make  the  arrival 
suite  difficult  to  interpret  accurately.  We  have  developed  a  set  of  techniques  for  small  aperture  arrays  that  improve 
phase  eharaeterization  and  identification  at  these  distances.  We  have  focused  our  efforts  on  two  arrays  in 
Kazakhstan  (MKAR  and  KKAR),  which  reeord  far-regional  events  throughout  South-Central  Asia.  Our  techniques 
include  improved  array  processing  methods  (e.g.,  phase- weight  semblance  stacking)  to  characterize  arrival  time, 
back-azimuth  and  slowness  of  individual  phases  within  the  P-coda  arrival  suite,  and  methods  to  more  accurately 
identify  the  measured  arrivals  (e.g.,  tau-p  transformation  and  clustering  analysis  to  determine  waveficld  templates). 

The  most  direct  means  ot  phase  identification  is  through  matching  arrival  time  and  slowness  estimates  to  theoretical 
predictions.  In  complex  tectonic  regions  such  as  South-Central  Asia,  global  reference  models  perform  poorly  at  far- 
regional  distances,  and  more  specific  know  ledge  of  the  along-path  and  near-array  earth  structure  is  required  for 
confident  primary  and  secondary  phase  identification.  To  account  for  near-array  effects  on  the  array  measurements, 
we  use  receiver  functions  to  image  below  the  arrays.  Receiver  functions  allow  a  direet  estimate  of  the  structure 
below  the  array  (e.g.  Moho  dip),  independent  of  the  array  measurements  and  without  bias  to  other  along-path 
effects.  We  then  use  array  measurements  (slowness  and  arrival  times)  to  empirically  derive  regionalized 
velocity-depth  profiles  that  more  accurately  predict  the  far-regional  phase  succession.  These  are  based  on  tau-p 
transformation  of  both  array  measurements  and  waveforms  to  form  generalized  representations  of  the  along-path 
earth  structure  that  is  sensitive  to  far-regional  propagation.  To  understand  the  far-regional  phase  behavior  not 
predicted  by  tau-p  ray  methods,  we  use  the  vcloeity-depth  profiles  to  synthesize  suites  of  seismograms  that  arc  used 
as  part  of  a  w  aveform  clustering  algorithm.  The  clustering  algorithm  processes  array  beams  to  derive  ’waveficld 
templates',  i.c.,  grouped  observations  with  similar  phase  characteristics.  The  waveficld  templates  are  further 
analyzed  in  a  non-linear  fashion  by  comparing  them  with  the  synthetic  seismograms,  looking  for  quantitative 
explanations  for  the  phase  behaviors  not  predicted  by  tau-p  ray  methods.  This  w  ork  has  resulted  in  a  methodology 
that  improves  the  usefulness  of  small-aperture  array  in  the  phase  eharaeterization  of  far-regional  earthquakes.  Our 
research  has  also  yielded  insight  into  body-wave  phases  that  arc  regularly  observed  on  the  MKAR  and  KKAR 
arrays,  including  information  on  expected  wave  propagation  behavior  and  the  regional  nature  of  the  upper-mantle 
discontinuities. 
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OBJECTIVE 

The  objective  of  our  research  is  to  enhance  the  usefulness  of  regional,  small-aperture  arrays  by  developing  array- 
based  methods  that  can  more  accurately  characterize  far-regional  (14  -29  )  P- coda  arrivals.  At  these  distances 
seismograms  are  particularly  interesting,  due  to  multipathing  through  the  lithosphere  and  upper-mantle.  A  suite  of 
body-wave  arrivals  is  often  observed  in  a  short  time  window  (1-20  seconds)  that  bottom  near  the  discontinuities  at 
220  km,  410  km,  and  660  km  depth,  depending  on  epieentral  distance  and  along-path  structure.  The  arrivals  can  be 
closely  spaced  in  time  and  exhibit  amplitude  anomalies  related  to  triplication  effects.  When  these  arrivals  can  be 
properly  identified  they  can  be  used  to  further  enhance  information  regarding  seismic  source  parameters,  as  well  as 
earth  structure.  However,  far-regional  seismograms  exhibit  significant  complexity  due  to  regional  heterogeneity  and 
interference  effects  from  depth-phase  multiples  and  near-reeeiver  scattered  arrivals.  It  is  the  goal  of  this  project  to 
improve  the  characterization  of  the  far-regional  arrivals  from  South-Central  Asian  events  by  employing  refined 
array  processing  techniques. 

The  regional  seismic  arrays  that  have  been  built  in  the  last  fifteen  years  should  be  a  rich  data  source  for  the  study  of 
far-regional  phase  behavior.  The  arrays  are  composed  of  high-quality  borehole  seismometers  that  make  high  fidelity, 
low-noise  recordings.  However,  beyond  regional  distances,  the  small  apertures  of  these  arrays  (<  5km)  limit  their 
usefulness  beyond  first-arrival  P-  and  S-  onset  picks.  Standard  array  methods  (e.g.,  slant  stacking  and 
frequency- wavenumber  analysis)  poorly  resolve  the  azimuth  and  slownesses  of  primary  and  secondary'  arrivals, 
making  confident  arrival  identification  and  classification  quite  difficult.  We  are  overcoming  this  limitation  by 
applying  refined  array  processing  techniques  in  conjunction  with  straightforward  wavefield  generalization 
methodologies.  Our  goal  is  to  characterize  the  commonly  observed  early  coda  arrivals  that  propagate  from  the 
different  seismic  regions  of  South-Central  Asia,  utilizing  recordings  from  the  Makanehi  (MKAR)  and  Karatau 
(KKAR)  arrays  in  Kazakhstan  (Fig.  1).  This  research  improves  the  usefulness  of  small-aperture  array  s  by  increasing 
their  ability  to  classify  small-magnitude  events  that  may  be  poorly  recorded  regionally  and  teleseismieally. 

RESEARCH  ACCOMPLISHED 

We  have  developed  array-based  methodologies  to  characterize  the  P  coda  of  far-regional  ev  ents  using  small  aperture 
arrays.  Our  methodology  includes  several  components:  1)  enhanced  array  calibration  to  account  for  near-array 
structure;  2)  improved  small-aperture  array  processing  to  effectively  measure  delay-times  (t)  and  slownesses  ip)  of 
primary  and  secondary  arrivals;  3)  construction  of  region-specific  velocity-depth  profiles  derived  from  the  t -p 
measurements  (i,e.,  x -/?  transformations),  as  well  as  models  derived  from  wavefield  continuation  methods;  and  4)  the 
derivation  of  ‘wavefield  templates’  from  clustering  of  array  beams  to  capture  the  commonly  observ  ed  arrival 
structure.  These  methods  employ  both  theoretical  and  empirical  techniques  to  reduce  a  large  waveform  dataset 
(-600  earthquakes;  Fig.  1 )  into  a  smaller  set  of  robustly  observed  wave  phenomena.  We  then  attempt  to  explain 
these  phenomena  using  well-accepted,  straightforward  techniques.  The  separate  components  of  our  methodologies 
are  discussed  below'. 
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Figure  1.  Map  of  South-Central  Asia  show  ing  the  location  of  the  Makanehi  (MKAR)  and  Karatau  (KKAR) 
arrays  in  Kazakhstan,  as  well  as  earthquakes  (circles)  used  in  this  study.  Colored  bands  mark  the 
14‘-29  distance  range  from  each  array. 
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Near-Array  Earth  Structure  from  Receiver  Functions 

Array-based  observations  of  slowness  and  back-azimuth  typically  deviate  from  ID  earth  predictions  because  of 
structural  heterogeneity  along  the  ray  path.  These  deviations  need  to  be  corrected  in  order  to  confidently  identify 
particular  arrivals.  In  many  cases,  heterogeneities  directly  below  the  array  can  account  for  most  of  the  observed 
anomalies  (c.g.  Niazi,  1966;  Havskov  and  Kansewich,  1978).  Interactions  with  upper  mantle  structure,  such  as 
topographic  relief  on  the  discontinuities,  can  also  produce  anomalous  observations.  At  the  MKAR  and  KKAR 
arrays,  we  observe  both  back-azimuth  and  slowness  anomalies,  where  azimuth  residuals  arc  ±20°  in  some  instances. 
We  think  that  near-array  structure,  as  well  as  along-path  effects,  contribute  to  these  observations. 

A  common  method  to  account  for  near-array  structure  is  to  calibrate  the  array  by  correcting  the  array-based 
measurements  to  theoretical  values  of  slowness  and  back-azimuth.  For  example,  published  earthquake  epicenters  arc 
used  to  determine  the  expected  great-circle  azimuth  angle,  and  a  ID  earth  model  is  used  to  compute  expected  ray 
parameters.  In  the  case  of  a  single  non-lateral  interface,  (e.g.,  dipping  Moho),  the  azimuth  and  slowness  residuals, 
relative  to  the  theoretical  values,  have  a  sinusoidal  pattern  with  respect  to  back-azimuth.  The  amplitude  and  phase  of 
these  patterns  can  be  used  to  solve  for  dip  and  strike  of  the  interface,  and  thus  a  correction  can  be  computed  for  each 
event  based  on  its  observ  ed  back-azimuth  angle.  This  method  works  well  only  if  near-array  structure  causes  the 
anomalous  residuals,  the  residuals  are  systematic  with  respect  to  back-azimuth,  and  there  is  good  azimuth  coverage 
of  array  measurements.  However,  for  the  MKAR  and  KKAR  arrays  these  criteria  are  not  well  met,  and  wc  require  a 
different  approach. 

To  obtain  a  more  accurate  accounting  of  the  ncar-array  structure,  and  thus  improve  our  array  measurements,  we  are 
employing  receiver  function  techniques  to  image  the  structure  below  the  arrays.  Our  goal  is  to  separate  near-array 
effects  on  the  azimuth  and  slowness  measurements  from  the  far-off  path  effects  by  accounting  for  the  ncar-array 
structure  independent  of  the  array  measurements.  Determining  structure  from  receiver  functions  has  several  benefits 
over  array-based  methods:  I )  only  the  ncar-array  structure  is  imaged  and  far-off  path  effects  have  no  influence,  2) 
earth  structure  is  not  determined  from  the  measurements  you  are  trying  to  correct  (the  array  measurements  in  our 
case),  and  3)  teleseismic  events  used  for  the  receiver  functions  arc  more  azimuthal ly  distributed  compared  to  the 
far-regional  events. 

For  the  MKAR  and  KKAR  arrays,  wc  compute  receiver  functions  using  data  from  the  single  3-component 
instrument  installed  at  each  array  using  a  time  domain  deconvolution  method.  The  distribution  of  teleseismic  events 
is  shown  in  Figure  2,  displayed  as  the  Moho  P-lo-S  conversion  points  around  each  array.  Several  hundred  receiver 
functions  arc  computed  for  each  array.  The  distribution  for  the  MKAR  array  is  more  azimuthally  complete  than  the 
KKAR  array,  but  wc  are  attempting  to  acquire  more  events  to  fill  the  gaps.  Receiver  function  images  arc  shown  in 
Figure  3  for  MKAR  and  Figure  4  for  KKAR. 


MKAR  KKAR 


Figure  2.  Azimuthal  distribution  of  teleseismic  events  used  in  the  receiver  function  calculations.  The  green 
dots  mark  the  Moho  conversion  point  for  each  of  the  receiver  functions. 
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Using  the  radial  and  tangential  receiver  functions  we  form  azimuthal  record  sections  to  image  P-Xo-S  converted 
phases  below  each  array.  The  record  sections  (Figures  3  and  4)  arc  constructed  by  stacking  receiver  functions  at 
constant  azimuth  spacing.  To  account  for  variation  in  ray  parameter  of  the  tclcscismic  events  and  to  minimize 
incoherent  scattered  energy,  we  use  a  2D  Gaussian  weighting  function  where  the  weight  is  a  function  of  nearness  to 
the  back-azimuth  stacking  point  and  ray  parameter  (Neal  and  Pavlis,  1999).  In  this  manner,  we  are  able  to  construct 
an  evenly  sampled  (in  back-azimuth),  smoothed  record  sections  that  can  be  used  in  subsequent  analysis. 

MKAR  Radial  MKAR  Tangential 


Figure  3.  MKAR  receiver  functions  record  sections.  The  receiver  functions  are  aligned  on  the  minimum- 
phase  time  and  have  been  filtered  at  0.1  to  0.4  Hz. 

Figures  3  and  4  show  the  MKAR  and  KKAR  record  sections.  The  main  P-\o-S  signal  at  both  MKAR  and  KKAR  has 
a  delay  of  ~6.0  seconds  for  the  high-frequency  band  we  arc  examining.  It  is  most  obvious  in  the  radial  component 
sections,  and  it  shows  some  variation  in  arrival-time  and  amplitude  for  different  back-azimuths.  Preceding  this  phase 
is  a  large  amplitude  negative  pulse,  which  may  be  related  to  a  low-velocity  zone  near  the  surface.  Later  arriving 
reverberations  are  also  apparent  in  both  the  MKAR  and  KKAR  sections  and  show  variations  with  back  azimuth 
Coherent  energy  is  apparent  on  the  tangential  component  sections,  which  is  typically  indicative  of  dipping  structure 
or,  in  some  cases,  anisotropy.  Coherent  energy  on  the  tangential  components  is  clearest  in  KKAR  record  sections, 
where  amplitude  polarity  changes  arc  also  seen  However,  observed  phases  and  the  coherency  of  the  record  sections 
seem  to  have  frequency  dependence.  At  lower  frequencies,  the  main  P-Xo-S  phase  is  less  prominent  and  the 
tangential  record  sections  arc  less  coherent  (not  shown).  This  may  be  related  to  the  nature  of  the  velocity  contrasts 
across  the  discontinuities  below  the  arrays,  and  will  likely  be  sorted  out  from  our  analysis 

We  are  still  in  the  process  of  analyzing  and  interpreting  the  receiver  function  images.  The  movcout  of  the  radial 
component  reverberations  and  the  existence  of  coherent  energy  on  the  tangential  components  receiver  functions 
seem  to  imply  that  the  earth  structure  below  the  arrays  is  not  laterally  homogeneous,  as  we  expected  from  the  array 
measurements.  In  order  to  understand  the  earth  structure  we  are  taking  a  simple  forward  modeling  approach.  We 
start  with  a  simple  flat-layer  velocity  model  determined  from  inverting  a  cumulative  stack  of  the  radial  receiver 
functions.  Using  this  model,  we  systematically  vary  the  strike  and  dip  of  the  Moho  to  compute  a  suite  of  synthetic 
receiver  functions.  The  synthetic  radial  and  tangential  record  sections  arc  then  quantitatively  compared  to  the 
observed  record  sections.  In  this  manner  we  arc  trying  to  solve  for  the  strike  and  dip  that  minimizes  the  difference 
between  the  observed  and  computed  record  sections.  Preliminary  results  from  this  exercise  for  the  MKAR  array 
show  that  much  of  the  receiver  function  arrival  structure  can  be  explained  by  a  prominent  discontinuity  at  48  km 
depth  striking  2 1 0°  N  with  a  dip  of  1 0°. 
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Figure  4.  KKAR  receiver  function  record  sections.  The  receiver  functions  arc  aligned  on  the  minimum-phase 
time  and  have  been  filtered  at  0.1  to  0.4  Hz. 


Improved  Small-Aperture  Array  Processing 

Our  objective  for  the  array  processing  is  to  compile  a  set  a  phase  delay  times  (x),  back-azimuth  and  slowness 
measurements  (/>)  and  array  beams  for  our  waveform  data  set  of -600  earthquakes.  The  small  aperture  (<  5  km)  of 
many  of  the  newer  arrays  installed  for  nuclear  monitoring  purposes  presents  a  serious  challenge  for  typical  array 
processing  methods  at  far-regional  distances.  The  restricted  array  aperture  and  limited  number  of  array  channels 
provides  limited  resolving  power,  particularly  for  slowness  measurements.  However,  we  have  developed  and 
applied  a  phase-weight  semblance  stacking  method  that  improves  the  usage  of  small-aperture  arrays  at  far-regional 
distances.  The  method  combines  an  amplitude  unbiased  measure  of  coherency,  the  phase-weight  stack  (e.g. 
Sehimmel  and  Paulssen,  1997),  with  semblance,  an  energy  eohereney  measure.  By  combining  the  two  coherency 
measures  the  signal-to-noise  level  is  increased  without  detrimentally  affecting  small  amplitude,  emergent  arrivals. 
The  semblance  stack  also  has  the  added  benefit  of  being  directly  related  to  the  F  statistic,  allowing  confidence 
testing  at  particular  signal-to-noise  levels.  Our  results  indicate  that  in  many  cases  we  can  detect  elosely-spaeed 
arrivals  by  their  slowness  values  in  a  time  window  that  includes  the  direct  P  arrival,  depth  phases  and  arrivals  from 
upper-mantle  discontinuities  (Ferris  and  Reiter,  2007). 

In  Figure  5  we  show  two  examples  of  applying  the  phase-weight  semblance  stacking  to  far-regional  events  recorded 
by  the  MKAR  array.  The  first  event  was  recorded  at  a  distance  of  25°  where  the  main  P  arrivals  bottom  near  the  660 
km  discontinuity  and  appear  in  the  first  few  seconds  of  the  seismograms.  The  other  dominant  phases  are  the/>P  and 
sP  depth  phase  packets  that  follow  a  similar  path,  exhibiting  similar  slowness  structure  as  the  main  P  arrivals  (as 
seen  in  the  slowness  vespagram)  but  increased  signal  amplitude.  The  second  event  was  recorded  at  a  distance  14.9°. 
At  this  distance  the  first  arrivals  are  Pn  followed  by  arrivals  bottoming  near  220  km  depth.  Later  arrivals  include  the 
P4 1 0  phase  (-8  seconds  later),  followed  by  larger  amplitude  7*41 0  depth  phases. 

The  back  azimuth  vespagrams  show  variations  in  azimuth  over  the  20-second  arrival  window  of  the  P  phases.  This 
is  most  ev  ident  for  the  event  at  25°  distance  and  reflects  out  of  plane  effects  from  along-path  and  near-array 
heterogeneity.  While  the  slowness  vespagrams  are  able  to  isolate  the  individual  arrivals  as  a  function  of  time,  the 
phase  weight  semblance  shows  smearing  in  the  slowness  dimension.  This  is  a  related  to  wave-number  aliasing  and  is 
a  function  of  the  array  configuration.  Further  improvement  to  the  slowness  and  azimuth  vespagrams  may  only  be 
realized  by  increasing  the  number  of  sensors  at  the  arrays  and  modifying  its  geometry. 
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Figure  5.  Example  of  phase-weight  semblance  processing  for  two  separate  earthquakes  recorded  by  the 

MKAR  array  at  a  cpiccntral  distance  of  25°  (left)  and  14.9°  (right).  The  top  row  shows  the  array 
gathers,  the  middle  row  shows  the  arrival  back-azimuth  vespagrams,  and  the  bottom  row  shows  the 
slowness  vespagrams. 


As  part  of  the  slowness  and  baek-azimuth  measurement  exercise  we  compared  our  slow  ness  estimates  to  those 
predicted  by  the  AK135  reference  model.  This  illuminated  the  need  for  more  region-specific  models  in  order  to 
make  accurate  identifications.  The  next  component  of  our  improved  phase  characterization  sehemc  uses  the  array 
measurement  and  subsequent  array  beams  to  derive  models  that  can  more  accurately  mateh  the  observations.  Our 
methodologies  to  accomplish  this  and  results  are  discussed  in  the  following  sections. 

Regional  Velocity  Depth  Profiles 

To  address  the  need  for  more  applicable  earth  models  and  to  generalize  our  array  observations,  we  derived 
veloeity-depth  profiles  for  the  specific  regions  monitored  by  the  MKAR  and  KKAR  arrays.  We  employed  two 
methods  based  on  r-p  (i.e.  delay-timc/slowncss)  techniques  that  decompose  wavefields  into  their  components  parts. 
The  first  method  involves  applying  wavcfield  continuation  techniques  (i.e.,  x-p  transformations)  to  the  array-based 
delay-time  and  slowness  measurements  to  solve  for  a  veloeity-depth  function  through  a  migration  procedure.  The 
second  method  applies  the  same  techniques  to  full-waveform  reeord  sections,  accomplishing  the  same  goal  (c.g., 
MeMcchan,  1984).  While  these  methods  are  not  new,  they  arc  not  typically  applied  to  small-aperture  array  data,  but 
rather  to  data  from  long  linear  arrays  of  seismometers  (e.g.,  Morozov  ct  al.,  2005).  Since  we  apply  these  methods  to 
multiple  events  records,  a  main  challenge  is  accounting  for  the  effects  of  differing  souree  parameters  (origin  time, 
event  depth,  and  focal  mechanism),  which  may  be  unknown  or  poorly  constrained.  Our  results  show  that  we  arc  able 
to  extract  reasonable  models  for  specific  regions  in  our  study  area  that  can  more  consistently  match  observed  arrival 
parameters. 

Figure  6  illustrates  how  we  prepare  array-based  delay-time  and  slowness  measurements  for  velocity-depth  inversion. 
In  this  example  we  groom  the  outliers  from  data  set  of  measurements  made  at  MKAR,  consisting  of  1 1 0  earthquakes 
ranging  from  14°-29°  epicentral  distance,  which  includes  earthquakes  extending  from  the  Hindu  Kush  region  of 
Pakistan/Afghanistan  to  the  Makran  coast  and  Zagros  Mountains  of  Iran.  The  grooming  removes  everything  outside 
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the  8.0  -  14.5  (sec/deg)  slowness  range,  to  help  eliminate  measurements  from  P-Xo-S  scattered  signals  and  coherent 
noise.  We  also  discard  measurements  that  exhibit  significant  slowness  smearing  (i.e.,  large  measurement 
uncertainties)  during  the  processing.  Once  the  groomed  measurements  (Figure  6a)  have  been  collected,  we  reduce 
them  to  a  single  r -p  curve  by  computing  the  error-weighted-mean  r  value  within  evenly  spaced  slowness  bins 
(Figure  6b),  excluding  bins  with  only  a  few  measurements.  Averaging  the  r’s  within  a  common  slowness  bin  acts  to 
correct  for  errors  in  focal  depth  by  accounting  for  the  ±r  offset  between  P  and  pP.  While  this  is  not  a  perfect 
correction  due  to  missed  phases  in  the  array  processing,  the  inclusion  of  sP  arrivals,  and  measurement  error,  we  have 
found  it  more  effective  and  feasible  than  correcting  r  on  an  event  basis  using  catalog  depths,  whieh  typically  have 
large  uncertainties.  The  averaged  r-p  curves,  while  useful  in  their  own  right,  are  further  processed  in  velocity-depth 
using  wavefield  continuation  methods,  which  we  describe  below'  as  applied  to  full-waveform  record  sections. 

The  wavefield  continuation  method  consists  of  two  linear  transformations  of  the  wavefield.  These  are  illustrated  in 
Figure  7.  First,  record  sections  are  slant  stacked,  transforming  them  from  the  distance-time  domain  into  the  r-p 
domain  (Figure  7b).  This  is  followed  by  a  downward  continuation,  or  depth  migration,  of  the  wavefield  to  transform 
the  r-p  data  to  the  slowness-depth  plane.  This  second  step  is  an  iterative  process,  where  the  i -p  data  are  repeatedly 
migrated  until  the  slowness-depth  image  converges  to  the  input  velocity  model  (Figure  7e).  There  are  several  ways 
to  update  the  velocity  model  between  iterations.  We  are  currently  using  a  scheme  that  picks  the  maximum  amplitude 
at  each  depth  from  the  slowness-depth  image  and  then  computing  the  weighted  average  between  it  in  the  input 
velocity  model.  This  scheme  performs  adequately;  however,  noise  in  the  initial  i-p  transformation  can  cause 
artifacts  in  the  final  velocity-depth  profile  that  needs  to  be  corrected  a  posteriori. 
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Figure  6.  Sample  r-p  data  from  the  MKAR  array,  a)  The  groomed  MKAR  data  (blue  circles);  red  dots  show 
the  iasp91  r-p  curve  for  a  surfaee-foeus  source.  The  range  in  r  for  a  particular  slowness  is  caused  by 
the  difference  in  r  for  earthquakes  at  different  depths,  rather  than  just  measurement  error,  b)  The 
averaged  r-p  eur\e  from  the  groomed  data,  evenly  sampled  in  slowness  at  0.15  see/deg. 


The  resulting  velocity-depth  profiles  show  general  agreement  between  the  array-based  measurement  transforms  and 
full-waveform  transforms.  The  models  more  accurately  match  the  far-regional  arrival  observations  at  MKAR  and 
KKAR  than  reference  models  such  as  iasp91  and  AKI35.  The  main  upper-mantle  discontinuities  are  found  at  250 
km,  415  km  and  670  km  depth  (Fig  7).  In  general,  the  wavefield  continuation  results  exhibit  moderately  slower 
velocities  than  iasp91  below  250  km  depth,  and  a  gradient  zone  is  observed  near  410  km  depth  rather  than  a  sharp 
discontinuity. 
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Figure  7.  Example  showing  the  transformation  of  a  record  section  (a)  to  the  r-p  domain  (b)  and  the 

downward  continuation  to  the  slowness-depth  domain  (c).  Array  beams  are  from  104  earthquakes 
that  extend  from  the  Hindu  Kush  region  of  Pakistan/Afghanistan  to  the  Makran  coast  and  Zagros 
mountains  of  Iran. 


Waveform  Templating  for  Groups  of  Like  Waveforms 

Another  component  of  our  phase  characterization  methodology  is  the  analysis  of  the  array  beams  using  clustering 
techniques.  In  addition  to  the  arrival  times  and  slowness  of  the  coda  phase  arrivals,  phase-weight  semblance 
stacking  analysis  yields  array  beams,  which  are  then  input  to  a  waveform  clustering  algorithm.  The  algorithm  is 
based  on  ‘fuzzy  clustering’  to  form  groups  of  beams  that  exhibit  similar  characteristics;  it  has  been  used  widely  in 
pattern  recognition  applications  (Bezdek,  1981 ).  Clusters  arc  defined  by  assigning  each  beam  a  cluster  membership 
value,  which  is  a  quantitative  measure  that  incorporates  a  distance  measure  between  each  array  beam  and  a 
representation  of  the  cluster  centers.  We  have  experimented  with  distance  measures  involving  LI  and  L2  norms,  as 
well  as  waveform  semblance,  with  varying  degrees  of  success  that  depend  on  data  noise  levels.  From  each  cluster  of 
similar  beams,  we  derive  a  template  waveform  using  a  stacking  process  that  weights  each  cluster  member  by  its 
degree  of  membership  to  that  particular  cluster.  The  ‘fuzziness’  of  the  method  allows  a  single  beam  to  be  a  member 
of  more  than  one  cluster  group.  The  objective  behind  the  waveform  clustering  is  to  reduce  the  database  of 
observations  to  a  set  of  representative  waveforms  that  exhibit  consistent  phase  arrival  behavior.  Sometimes  the 
clustering  is  geographic,  but  there  are  other  wave  phenomena  that  can  also  produce  groupings  of  similar  waveforms 

Figure  8a  shows  the  cluster  group  and  its  associated  wavefield  template  that  results  from  applying  our  waveform 
clustering  algorithm  to  the  set  of  100  events  shown  in  Figure  8b.  These  beams  arc  derived  from  seismograms 
recorded  by  the  MKAR  array  in  the  14°-29°  distance  range  that  extends  from  northern  Pakistan  to  the  Makran  coast 
in  Iran.  In  this  example  the  algorithm  clusters  the  events  into  5  groups,  with  12-24  members  per  group.  We  applied 
an  Z,2-norm  distance  measure  and  aligned  the  array  beams  prior  to  clustering  using  multi-channel  cross  correlation 
applied  to  the  first  four  seconds  of  each  beam  signal.  We  have  found  the  alignment  of  the  beams  to  be  beneficial  to 
the  clustering  procedure.  However,  the  length  of  the  eross-eorrelation  window  is  an  important  variable.  For  example, 
if  the  window  is  too  long,  the  array  beams  align  on  the  maximum  amplitude  signal  in  the  record.  For  some  beams 
the  maximum  amplitude  arrival  may  occur  early  in  the  seismogram,  while  for  others  it  occurs  much  later.  This  can 
result  in  some  unusual  clustering  of  events  that  seems  counter-intuitive.  For  this  reason,  we  find  clustering  methods 
based  strictly  on  cross  correlation,  as  used  in  other  types  of  studies  (e.g.,  Menke  1999)  are  not  suitable  for  our 
particular  needs. 
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Figure  8.  Example  of  wavefield  clustering  to  generate  a  wavefield  template,  (a)  The  bottom  panel  shows 
the  cluster  members  sorted  by  epiccntral  distance  (shown  on  the  left).  The  green  band  highlights 
the  time  band  of  cross  correlation.  The  top  panel  shows  the  computed  w  avefield  template  (black) 
and  1-<t  deviation  (red),  (b)  A  record  section  of  all  100  earthquakes  input  to  the  clustering 
algorithm. 


For  the  example  shown  in  Figure  8,  the  resulting  cluster  members  include  earthquakes  that  are  near  the  14.3° 
distance  range,  but  some  more  distant  earthquakes  are  also  included.  The  earthquake  depth  range  is  also  variable, 
but  may  be  due  to  poorly  constrained  catalog  depths.  In  general,  the  cluster  members  show  consistent  waveform 
structure  in  the  cross-correlated  portion  of  the  signal.  Later  arrivals  show  more  variability  between  the  cluster 
members,  but  still  appear  to  align  well  at  approximately  8  seconds  and  later.  The  most  variability  occurs  between 
4  and  8  seconds  delay  time,  where  some  records  show'  large  amplitude  arrivals  while  other  show  little  signal.  Some 
of  this  variability  is  likely  due  to  differences  in  earthquake  distance,  triplicated  arrivals,  and  the  presence  of  depth 
phases.  We  continue  to  refine  and  improve  the  wavefield  clustering  algorithm. 


Phase  Characterization  Analysis 

A  final  component  of  our  far-regional  methodology  is  to  characterize  the  observed  phase  arrivals.  One  way  we  have 
been  attempting  to  do  this  is  to  quantitatively  compare  arrivals  predicted  from  the  derived  velocity-depth  profiles  to 
those  observed  in  the  ‘wavefield  templates’.  The  purpose  of  this  exercise  is  to  determine  how  much  of  the  arrival 
structure  can  be  explained  by  simple  and  straightforward  methods  and  models.  Using  these  velocity-depth  models 
we  generate  a  suite  of  synthetic  seismograms,  which  w  e  then  use  to  compute  several  different  misfit  measures 
between  the  synthetics  and  the  observed  wavefield  templates.  Since  the  synthetic  waveforms  are  generated  from 
models  derived  from  the  x -p  data,  arrival  times  of  the  w  avefield  templates  should  be  w  ell  matched.  We  note  that  to 
derive  meaningful  measurements  we  must  account  for  frequency  content  and  amplitude  variations. 

While  this  exercise  in  waveform  modeling  is  relatively  crude,  we  are  finding  general  agreement  between  the 
synthetic  waveforms  and  the  ‘wavefield’  templates.  This  suggests  to  us  that  a  significant  portion  of  the  observed 
far-regional  arrival  structure  can  be  explained  by  regional-specific  ID  models. 

CONCLUSIONS  AND  RECOMMENDATIONS 

Regional  seismic  arrays  that  have  been  recently  installed  for  nuclear  monitoring  are  under-utilized  in  the  study  of 
far-regional  arriv  als.  The  small  aperture  of  many  of  these  arrays  (<  5km)  restricts  their  usefulness  at  these  distances 
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beyond  first  arrival  onset  picks  of  P  and  S  waves.  However,  our  research  is  overcoming  this  limitation  by  applying 
refined  array  processing  techniques  in  conjunction  with  x-p  and  wavefield  templating  analysis  methods.  Our 
methodology  improves  the  characterization  of  primary  and  early  coda  phase  arrivals  observ  ed  at  far-regional  and 
near-teleseismie  distances.  Our  approach  is  to  distill  the  wide  v  ariety  of  seismicity  we  observe  to  subsets  of 
commonly  and  robustly  observed  phase  arrival  phenomena.  We  are  then  using  well-aeeepted  modeling  and 
inversion  techniques  to  explain  these  phenomena.  Our  aim  is  to  explain  as  much  of  the  phenomena  as  we  can  with 
simple  and  straightforward  techniques,  leav  ing  the  anomalies  for  future  research 

We  are  developing  and  applying  our  techniques  to  South-Central  Asian  earthquakes  recorded  on  the  MKAR  and 
KKAR  arrays  in  Kazahkstan.  Our  results  indicate  that  we  can  differentiate  between  the  numerous  arrivals  of  the 
early  P-eoda.  However,  global  reference  models  cannot  capture  the  phase  succession  and  arrival-time  behavior  we 
observe  from  the  complex  tectonic  regions  of  South-Central  Asia  To  address  this,  we  have  developed  regional  ID 
models  directly  from  the  array  measurements  (delay-times,  slowness,  and  beams).  Since  these  models  are  derived 
from  the  data,  they  are  able  to  explain  the  phase  behavior  we  observe  from  specific  regions.  To  test  these  models 
and  further  improve  phase  characterization,  we  have  constructed  ‘wavefield  templates’  through  cluster  analysis  to 
generalize  the  waveform  structure  from  the  difference  seismic  regions.  The  templates  are  then  used  in  a  waveform 
fitting  analysis  to  gain  a  better  understanding  of  the  phase  phenomenon  observed  from  the  complex  seismic  regions 
of  central  Asia. 

While  our  methodology  improves  the  usefulness  of  small-aperture  arrays  at  analyzing  far-regional  P-eoda  arriv  als, 
greater  gains  could  be  made  by  installing  more  elements  and  increasing  the  effective  array  aperture  at  both  MKAR 
and  KKAR.  Modifying  the  array  geometry  would  reduce  the  wave-number  aliasing  for  arrivals  with  slowness  values 
between  8-1 1  s/deg  and  increasing  the  aperture  would  improve  both  slowness  and  azimuth  resolution.  An  addition  of 
2-6  new  sensors  in  an  outer-ring  configuration  would  significantly  improve  far-regional  array  analysis. 
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